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ABSTRACT 


This  report  represents  the  second  chapter  of  a  book  in  preparation  on 
inference  and  data  analysis  in  reliability  and  life  testing.  The  point 
of  view  adopted  differs  from  that  of  most  books  on  the  subject  in  the 
following  basic  respect:  Prior  information  available  to  the  reliability 
analyst  is  utilized  fully  in  a  formal  statistical  fashion.  Experience 
accumulated  in  helping  engineers,  quality  assurance  managers,  scientists, 
biostatisticians,  and  others  who  must  make  estimates  and  reach  decisions 
from  either  planned  experiments  or  retrospective  data  has  shown  us  that 
the  point  of  view  adopted  throughout  the  book  has  resulted  in  useful 
solutions  to  real-life  problems.  By  contrast,  more  classical  statistical 
methods  have  often  proven  inadequate  in  many  practical  problems  simply 
because  the  data  available  are  insufficient  to  reach  conclusions  with  a 
desired  degree  of  assurance. 

The  book  is  intended  primarily  for  actual  use  by  the  engineering  and 
scientific  practitioner,  rather  than  for  theoretical  study  and  philo¬ 
sophical  analysis  by  the  statistician.  Thus  we  omit  a  philosophical 
justification  of  the  methods  presented;  rather,  we  rely  on  the  fact  that 
they  have  led  to  useful  answers  to  problems  that  have  arisen  in  practice. 


One  final  point:  Many  of  the  methods  and  results  are  original  and  have 
not  appeared  in  the  literature.  This  fact  has  led  us  to  issue  the 
chapters  as  reports  under  our  research  grants. 
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LIFE  DISTRIBUTION  MODELS  AND  INCOMPLETE  DATA 


by 

Richard  E.  Barlow  and  Frank  Proschan 

In  this  chapter  our  objective  is  to  introduce  additional  life  dis¬ 
tribution  models  and  to  discuss  methods  useful  for  analyzing  failure 
data,  especially  incomplete  data.  We  show  how  to  express  the  like¬ 
lihood  functions  for  general  distributions  and  incomplete  data.  The 
likelihood  function  tends  to  be  fairly  flat  for  incomplete  data.  For 
this  reason  the  maximum  likelihood  estimator  may  be  of  limited  value. 

It  is  therefore  especially  important  in  this  situation  to  assess  a 
prior  distribution  for  parameters  and  plot  the  posterior  distribution 
or  its  contours. 

Inference  based  on  the  exponential  model  is  discussed  for  general 
sampling  plans.  Parameter  estimators  and  credibility  intervals  are 
derived  for  special  cases.  The  Weibull  distribution  is  a  very  useful 
model  for  life  distribution  studies  and  also  for  the  analysis  of 
strength  data.  For  these  reasons  we  describe  failure  mechanisms 
leading  to  a  Weibull  life  distribution  model.  Contour  plotting  methods 
for  analyzing  life  data  based  on  a  Weibull  distribution  are  also  given. 

1.  LIKELIHOOD 

In  Chapter  1,  we  calculated  the  likelihood  for  the  exponential 
model  under  several  different  sampling  plans.  In  this  section  we 
present  a  unified  way  of  analyzing  incomplete  data  for  a  large 
number  of  failure  distribution  models.  In  much  of  this  chapter  we 
assume  that  the  failure  distribution  F  is  absolutely  continuous 
with  density  f  and  failure  rate 
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(1.1) 


r(x) 


f(x) 
F(x)  ’ 


where  F(x)  =  1  -  F(x)  .  We  call 


X 

(1.2)  R(x)  =*  j  r(u)du 

0 

the  hazard  function  associated  with  F  .  For  general  F  ,  define 

(1.3)  R(x)  *  -In  F(x) 

so  that  F(x)  *  exp  [-R(x) ]  .  Note  that  when  F  has  a  density  f  , 

4:  [~ln  F(x)  ]  =  =  r(x) 

***  F(x) 

so  that  (1.2)  and  (1.3)  agree  in  this  case. 

From  (1.1)  and  (1.3)  we  see  that 

(1.4)  f (x)  -  r(x)e“R(x)  . 

For  a  discussion  of  these  fundamental  concepts,  their  inter-relation¬ 
ships  and  illustrations  in  the  case  of  well  known  distributions,  see 
Barlow  and  Proschan,  Chapter  3,  (1975). 

Suppose  now  we  observe  n  independent  lifetimes  x^,.:^,  ...»  xn 
corresponding  to  a  given  failure  rate  function,  r  .  The  joint  density 
is 


n  f(x.) 
i-1 


n 

n 

i-l 


r(xi) 


n 

-  I  R(x.) 

i-l 


(1.5) 
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The  likelihood  as  a  function  of  the  failure  rate  function  for  data 
D  *  . ..,  xn)  is  then 


(1.6) 


L(r (u)  ,  u  >  0  |  D)  = 


n 

n 

i-l 


r(x1) 


-  I  R(xt) 

i=l 


1.1  Example:  The  Time-Transformed  Exponential  Model 
Suppose  the  survival  function  is  of  the  form: 


-XR  (x) 

(1.7)  F(x  |  A)  -  e 

where  it  is  assumed  that  R  is  known  and  differentiable  but  X  is 

o 

unknown.  By  (1.2)  we  may  write: 

x 

XR  (x)  *  I  Xr  (u)du  . 

O  J  0 

0 

It  follows  that  the  hazard  function  and  the  failure  rate  function  are 

assumed  known  up  to  the  parameter  X  .  Another  way  to  view  the  model  is 

to  consider  time  x  to  be  transformed  by  the  function  R  (•)  .  Thus  in 

0 

terms  of  transformed  time,  the  present  model  coincides  with  the  exponen¬ 
tial  model  of  Chapter  1.  For  this  reason  (1.7)  is  called  the  time- 
transformed  exponential  model. 

Let  x.,x_,  ...,  x  be  n  independent  observations  given  X 
l  l  n 

from  this  model.  The  likelihood  is 


n 

n 

i-i 


ro(xi> 


"x  l ,R0(xi) 
1=1 


(1.8) 


L  (X  i  D)  =  An 


e 
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n 

By  Lemma  1.5  of  Chapter  1,  we  conclude  that  £  R  (x  )  and  n  is 

i-1  ° 

a  sufficient  statistics  for  X  .  If  we  use  the  gamma  prior  for  X  : 


ir(X) 


,a,  a-1  -bX 
b  X  e 

r(a) 


we  obtain  as  the  posterior  density  for  X  : 


(1.9)  tt(X  |  D) 


(See  (1.6)  of  Chapter  1.)  Inference  proceeds  exactly  as  for  the 
exponential  model,  except  that  observation  x^  of  the  exponential 
model  is  replaced  by  its  time-transformed  value  R^x^)  .  This  is 
valid  assuming  only  that  Rq(«)  is  continuous. 


The  General  Sampling 


In  many  practical  life  testing  situations,  the  lifetime  data 
collected  are  incomplete .  This  may  be  due  to  the  sampling  plan  itself 
(recall  the  sampling  plans  (a),  (b)  and  (c)  of  Chapter  1),  or  due  to 
the  unplanned  withdrawal  of  test  units  during  the  test.  (For  example, 
in  a  medical  experiment,  one  or  more  of  the  subjects  may  leave  town, 
or  suffer  an  accident,  etc.) 

We  now  describe  one  type  of  sampling  plan.  Suppose  unit  i 
having  lifetime  distribution  F  is  observed  over  an  interval  of  time 
starting  at  age  0  and  ending  at  a  random  or  nonrandom  age.  Termina¬ 
tion  of  observation  occurs  in  either  one  of  the  following  two  ways: 
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(1)  The  ith  unit  is  withdrawn  or  lost  from  observation  at  age 

>  0  ;  may  be  random  or  nonrandom. 

(2)  The  ith  unit  fails  at  age  X^  ,  where  is  a  random 

variable. 

In  addition,  we  require  a  technical  assumption  regarding  the 
"stopping  rule";  i.e.,  a  prescription  for  determining  when  to  stop 
observation: 

(3)  Suppose  unit  lifetime,  X  ,  depends  on  an  unknown  parameter 
(or  parameters)  9  .  Observation  on  a  unit  may  stop  before 
unit  lifetime  is  observed.  Let  S  be  a  rule  or  set  of 
instructions  which  determines  when  observation  of  a  unit  stops. 
S  may  be  random  and  depend  on  an  unknown  parameter  (or  para¬ 
meters)  <)>  .  Let  (9,<j>)  have  a  joint  prior  distribution. 
Suppose  that  the  rule  S  and  9  ,  given  <p  ,  are  judged  in¬ 
dependent.  Also  suppose  that  6  and  $  are  judged  indepen¬ 
dent. 

It  is  important  to  remark  that  the  "stopping  rule"  is  not  necessarily 
the  same  as  the  "stopping  time." 

To  understand  assumption  (3),  consider  sampling  plan  (a)  of 
Chapter  1,  namely:  put  n  items  on  life  test  and  stop  testing  at  the 
kth  observed  failure.  In  this  case,  the  stopping  rule  depends  only  on 
k  and  is  clearly  independent  of  life  distribution  parameters  since  k 
is  fixed  in  advance  of  testing. 

Sampling  plan  (b)  of  Chapter  1  has  stopping  rule:  stop  testing  at 
time  tQ  .  Since  tQ  is  fixed  in  advance  of  testing,  the  stopping 
rule  is  again  independent  of  life  distribution  parameters.  Likewise, 
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the  stopping  rule  for  sampling  plan  (c)  depends  only  on  k  and  t 

and  is  clearly  also  independent  of  life  distribution  parameters. 

For  all  three  sampling  plans,  the  likelihood,  up  to  a  constant  of 

proportionality,  depends  only  on  the  life  distribution  model  and  the 

observed  data.  This  proportionality  constant  depends  on  the  stopping 
n  ! 

rule  and  is  ^  — — jy  for  sampling  plans  (a)  and  (b)  while  for  sampling 

plan  (c)  it  is  ^  -^k)i(k^~  1)~1  '  Such  stopping  rules  are  said  to  be 
nonin  formative  since  the  posterior  distribution  in  each  case  depends 
only  on  the  life  distribution  model,  the  observed  data,  the  prior,  and 
in  these  examples,  on  the  independence  of  unit  lifetimes  given  distri¬ 
bution  parameters. 

1.3  Examples  of  Informative  Stopping  Rules 

Records  are  routinely  kept  on  failures  (partial  or  otherwise)  and 
maintenance  actions  on  critical  units  such  as  airplane  engines.  Should 
a  relatively  new  type  of  unit  start  exhibiting  problems  earlier  than  anti¬ 
cipated,  this  may  trigger  early  withdrawal  of  units.  If  this  happens 
the  stopping  rule,  which  is  contingent  on  performance,  may  also  be 
informative  relative  to  life  distribution  parameters.  This  fact  needs 
to  be  considered  when  calculating  the  likelihood  and  analyzing  the  data. 

The  second  example  illustrates  another  case  where  assumption  (3) 
is  violated.  Suppose  lifetime  X  is  exponential  with  failure  rate  X 
and  the  random  withdrawal  time,  W  ,  is  also  exponential  with  parameter 
<J>  .  We  observe  the  minimum  of  X  and  W  .  Furthermore,  suppose  that  X 
given  X  and  W  given  <j>  are  judged  independent.  Then  the  likelihood 
given  an  observed  failure  at  x  is 

L(X  ,  |  x)  -  Xe  ^Xe  ^X  . 
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If  X  and  <j>  are  judged  a  priori  independent  then  the  posterior  density 
of  X  is 

tt(X  |  x)  «  Xe  ^Xtt(X) 

where  it  is  the  prior  density  for  X  .  However,  if  X  and  $  are 
judged  dependent  with  joint  prior  ir(X,<j>)  ,  then  the  posterior  density  is 

00 

it (X  |  x)  «  Xe  J*  e  ^Xit(X  ,<j>)  d4>  . 

0 

00 

The  factor  J*  e  ^Xir(X ,<j>)d<ji  ,  contributed  by  the  stopping  rule,  depends 
0 

on  X  . 

There  is  an  important  case  not  covered  by  the  General  Sampling  Plan — 
namely  when  it  is  known  that  a  unit  has  failed  within  some  time  interval 
but  the  exact  time  of  failure  is  unknown.  This  case  will  be  covered  when 
we  discuss  life  tables. 

The  following  simple  example  illustrates  the  way  in  which  incomplete 
data  can  arise. 

1.4  Example 

Operating  data  are  collected  on  an  airplane  part  for  a  fleet  of  air¬ 
planes.  A  typical  age  history  for  several  engines  is  shown  in  Figure  1.1. 
The  crosses  indicate  the  observed  ages  at  failure.  Ordered  withdrawal 
times  (nonfailure  times)  are  indicated  by  short  vertical  lines.  In  our 
example,  units  2  and  4  fail  at  respective  times  x^  and  x(2)  while 
observation  on  units  1  and  3  is  terminated  without  failure  at  times  ^-(2) 
and  respectively. 
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Age  u 


X(l)  *(1)  X(2)  l(2) 

FIGURE  1.1 

AGE  OF  AIRPLANE  PART  AT  FAILURE  OR  WITHDRAWAL 


9 


It  is  important  to  note  that  all  data  are  plotted  against  the  age  axis 
Figure  1.2  illustrates  how  events  may  have  occurred  in  calendar  time.  For 
example,  units  1  and  3  had  not  failed  at  the  end  of  the  calendar  record. 

Total  Time  on  Test 

As  we  saw  in  Chapter  1,  the  total  time  on  test  is  an  important  statis¬ 
tic  for  the  exponential  model. 

1.5  Definition 

The  total  time  on  test  T  is  the  total  of  the  periods  of  observation 
of  all  the  units  undergoing  test.  Excluded  from  this  statistic  are  any 
periods  following  death  or  withdrawal  or  preceding  observation.  Specifi¬ 
cally,  the  periods  being  totalled  include  only  those  in  which  a  death  or  a 
withdrawal  of  a  unit  under  observation  can  be  observed. 

Let  n(u)  be  the  number  of  units  observed  to  be  operating  at  age  u  . 
The  observed  function  n(u)  ,  u  _>  0  ,  for  Example  1.4  is  displayed  in 
Figure  1.3.  From  Figure  1.3  we  may  readily  calculate  the  total  time  on 
test  T(t)  corresponding  to  any  age  t  ,  0  <  t  <  '• 

t 

(1.10)  T(t)  =  J"  n(u)du  . 

0 

For  example,  for  t  such  that  x(2)  <  c  <  "^(2)  ’  we  °^ta^n  from  Figure  1.3 
t 

T(t)  -  J*  n(u)du  -  +  3(£(i)  -  x(i))  +  2^x(2)  “  ^(1)^  +  ^  “  x(2)^  ‘ 

0 

After  simplifying  algebraically,  we  obtain: 


Start  of 
Calendar  Record 


End  of 

Calendar  Record 


FIGURE  1.2 


CALENDAR  RECORD  FOR  AIRPLANE  PARTS 


IUMBER  OF  UNITS  IN 
OPERATION  n(u) 


4  - 
3  - 
2  - 
1  - 

0 


Age  u 

X(2)  Z{2) 


FIGURE  1.3 

NUMBER  OF  UNITS  IN  OPERATION  AS  A  FUNCTION  OF  AGE 
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(1.11)  T(t)  “  x(l)  ^(1)  x(2)  ^  ' 

Note  that  the  resulting  expression,  given  in  (1.11),  can  be  obtained 
directly,  since  x^  and  x(2)  represent  the  observed  lifetimes  of 
the  2  units  that  are  observed  to  fail,  represents  the  observed 

age  of  withdrawal  of  the  unit  first  withdrawn  from  observation,  and 
finally  t  represents  the  age  of  the  second  unit  at  the  instant  t 
specified. 

Although  in  this  small  example,  the  directly  calculated  expression 
(1.11)  for  total  time  on  test  is  simpler,  Equation  (1.10)  is  an  important 
identity,  since  it  yields  the  total  time  on  test  accumulated  by  age  t 
in  terms  of  the  (varying)  number  of  units  on  test  at  each  instant  during 
the  interval  [0,t]  for  any  data  set  in  which  the  agea  at  death  or  with¬ 
drawal  are  observed.  Thus  it  is  a  general  formula  applicable  in  a  great 
variety  of  problems  in  which  data  may  be  incomplete. 

Although  n(u)  is  a  step  function,  the  integral  representation  in 
(1.10)  is  advantageous,  since  it  is  compact,  mathematically  tractable, 

and  applicable  in  a  great  variety  of  incomplete  data  situations.  Of 
00 

course,  J*  n(u)du  <  00  in  practical  problems  since  observation  ultimately 
0 

ceases  in  order  to  analyze  the  data  in  hand. 

The  Likelihood  Function  for  Incomplete  Data 

All  recorded  data  are  necessarily  discrete.  Likewise  real  world  life 
distribution  models  should  also  be  discrete.  Continuous  life  distribution 
models  are  convenient  approximations  to  real  world  life  distributions.  How¬ 
ever,  it  is  most  convenient  to  define  initially  the  likelihood  concept  in 


the  context  of  discrete  models. 
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For  our  purposes,  we  find  it  preferable  to  define  the  likelihood 
concept  for  the  General  Sampling  Plan  in  the  context  of  a  discrete 
model.  Computation  of  the  likelihood  function  is  an  intermediate  step 
between  specification  of  the  prior  distribution  on  the  space  0  and 
computation  of  the  posterior  distribution  on  0  given  observed  data  D  . 

Suppose  temporarily  that  the  life  distribution  is  discrete,  i.e., 
failures  can  occur  only  at  times  1,2,  ...  ;  similarly,  withdrawals  can 
occur  only  at  these  time  points.  Suppose  that  the  probability  of  failure 
of  a  given  unit  at  x  is  p(x  |  9)  .  Suppose  k  failures  are  observed 
at  times  x  ,  s  ■  1,  ...,k,  and  m  withdrawals  are  observed  at  times 
£  ,  t  -  1,  ....  m  .  Failure  and  withdrawal  times  need  not  be  distinct. 

All  observations  are  assumed  statistically  independent,  given  parameters. 
Withdrawal  times  are  produced  by  a  stopping  rule  which  is  noninformative 
concerning  0  . 

For  example,  the  stopping  rule  might  specify  that  we  observe  a  unit 
until  failure  or  until  withdrawal,  whichever  comes  first,  where  withdrawal 
time  is  specified  in  advance.  For  this  model,  the  probability  of  the 
observed  outcome  is: 

k  m 

(1.12)  P(D  |  e)  -  n  P(x  |  e)  n  PCL  |  9)  , 

s-l  t=l 

co 

where  P(u  |  9)  ■  £  p(u.  |  9)  represents  the  probability  that  a 

J  i-1  J  1 

specified  unit  fails  at  age  u  ^  or  later,  given  the  parameter  is 

0  .  Note  that  the  first  product  corresponds  to  the  k  failures  at 

respective  ages  x^,  ...,  x^  ,  while  the  second  product  corresponds 

to  the  m  withdrawals  at  respective  ages  ,  . . . ,  t 

1  m 

Another  way  to  model  withdrawal  is  to  suppose  there  exists  a 
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random  withdrawal  age  W  such  that  P[W  -  t]  *  q(t)  ,  t  *  1,2,  ...  , 
with  W  independent  of  unit  lifetimes  and  of  8  .  Under  this  model, 
we  suppose  that  we  observe 


minimum  (X,W) 


Now  for  observed  data  D  *  {x^,  ...,  x^  ,  t ...,  ,  the  probability 

of  the  observed  outcome  given  parameter  9  is: 


!x 

if 

X 

<_  w 

(w 

if 

X 

>  w  . 

/ 

*  *  * 

•  •  ’  A 

(1.13)  p(d  I  8)  -  n  q td  )  n  Q(x  )  n  P(x  |  9)  n  m,  I  0)  , 

t-1  s-1  s-1  t-1  c 


where  Q(u.)  £  q(u,  )  represents  the  probability  that  W  >  u.  . 

3  i»l  3  3 

Note  that  (1.12)  and  (1.13)  differ  only  by  a  factor  that  does  not  depend 
on  9  .  Thus,  relative  to  calculating  the  MLE  of  0  ,  the  two  models 
for  withdrawals  (withdrawal  deterministic  or  withdrawal  random)  do  not 
differ  essentially . 

There  are  many  practical  testing  situations  in  which  withdrawals 
occur  as  a  result  of  chance  mechanisms  unrelated  to  the  parameter  8 
of  the  lifetime  distribution.  For  example,  concluding  the  collection 
of  data  at  a  specified  chronological  time  has  the  effect  of  withdrawing 
from  observation  those  units  still  alive  at  that  point  in  time.  In 
Figure  1.2,  this  phenomenon  is  illustrated  by  units  1  and  3.  Other  chance 
mechanisms  causing  withdrawal  at  a  random  age  result  from  human  errors 
and  accidents.  The  net  effect  of  the  various  stopping  rules  that  are 
unrelated  to  the  value  of  the  parameter  8  is  summarized  in  the  factor 
g(x,£)  in  the  expression  for  the  probability  of  the  observed  outcome: 


(1.14) 


k  m 

p(D  I  9)  ■  g(x,£)  n  P (x  !  0)  n  ?U  1  e) 
S-l  J  t-1  c 
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1.6  Definition 

The  likelihood,  L(0  |  D)  ,  is  the  probability  of  the  observed 
outcome,  p(D  |  8)  ,  considered  as  a  function  of  the  parameter  9 
given  the  data,  D  .  In  the  case  of  a  continuous  model,  the  correspond¬ 
ing  likelihood  will  have  this  interpretation  relative  to  a  discrete 
probability  approximation. 

It  follows  from  (1.14)  that 

k  m 

(1.15)  L(0  |  D)  «  n  p(x  j  0)  n  P(l  |  8)  . 

s-1  t-1 

From  Bayes'  Theorem  (Chapter  1,  Theorem  1.2.),  it  is  clear  that  we  need 
not  know  g(x,£)  in  order  to  compute  the  posterior  density  of  9  . 

In  this  subsection,  we  have  thus  far  confined  our  discussion  to 
the  case  of  discrete  time  life  distributions  since  the  basic  concepts 
are  easier  to  grasp  in  this  case.  However,  in  the  case  of  continuous 
time  life  distributions,  the  likelihood  concept  is  equally  relevant, 
and  in  fact  the  expression  for  the  likelihood  L(0  j  D)  assumes  a 
rather  elegant  form  if  we  use  n(u)  ,  the  number  on  test  function. 

In  the  continuous  case,  p(x  )  0)  is  replaced  by  the  probability 
density  element  f(x  |  0)  . 

1.7  Theorem 

Given  the  failure  rate,  independent  observations  are  made  under  the 
General  Sampling  Plan.  Let  x^^,  •••,  denote  the  k  observed 
failure  ages.  Let  n(u)  denote  the  number  of  units  under  observation  at 
age  u  ,  u  >_  0  ,  and  r(u)  denote  the  failure  rate  function  of  the  unit 
at  age  u  .  Then  the  likelihood  of  the  failure  rate  function  r(u)  , 


L 


having  observed  the  data  D  described  above,  is  given  by: 


(1.16) 


L(r(u)  ,  u  >  0  D) 


r 

II  r(x  )  exp  -  I  n(u)r(u)du  ,  k  >_  1 
■-1  S  0 


wo 

-  j*  n(u)r(u 


)du  , 


k  =  0 


Proof : 


To  justify  (1.16),  we  first  note  that  the  underlying  random  events 
are  the  ages  at  failure  or  withdrawal.  Thus  the  likelihood  of  the 
observed  outcome  is  specified  by  the  likelihood  of  the  failure  ages 
and  survivals  until  withdrawal.  By  Assumption  (3)  of  the  General 
Sampling  Model,  we  need  not  include  any  factor  contributed  by  the  stop¬ 
ping  rule,  since  the  stopping  rule  does  not  depend  on  the  failure  rate 
function  r(*)  . 

To  calculate  the  likelihood,  we  use  the  fact  that  given  r(*)  , 


f(x)  -  r(x)e 


x 

-/ r(u)du 
0 


(See  (1.4).)  Specifically,  if  a  unit  is  observed  from  age  0  until 
it  is  withdrawn  at  age  without  having  failed  during  the  interval 


[0,£t]  ,  a  factor  exp 


'  *t 

-  J  r(u)du  i 
0 


is  contributed  to  the  like¬ 


lihood.  Thus,  if  no  units  fail  during  the  test  (i.e.,  k  ■  0)  ,  the 


likelihood  of  the  observed  outcome  is  proportional  to  the  expression 
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given  in  (1.16)  for  k  *■  0  . 

On  the  other  hand,  if  a  unit  is  observed  from  age  0  until  it  fails 
at  age  xg  ,  a  factor, 


r(xg)  exp 


5 

J 


r (u)du 


is  contributed  to  the  likelihood.  The  exponential  factor  corresponds  to 
the  survival  of  the  unit  during  [0,x  )  ,  while  r(x  )  represents  the 
rate  of  failure  at  age  xg  .  (Note  that  if  we  had  retained  the  differ¬ 
ential  element  "dx"  ,  the  corresponding  expression  r(xg)dx  would 
approximate  an  actual  probability:  the  conditional  probability  of  a 

failure  during  the  interval  (x  ,  x  +  dx)  given  survival  to  age  x  .  ) 

s  s  s 

The  likelihood  expression  in  (1.16)  corresponding  to  the  outcome 
k  >_  1  now  is  clear.  The  exponential  factor  corresponds  to  the  survival 
intervals  of  both  units  that  failed  under  observation  and  units  that 
were  withdrawn  before  failing: 


J 


n(u)r (u)du 


r(u)du  + 


r (u)du 


where  the  first  sum  is  taken  over  units  that  failed  while  the  second  sum 
is  taken  over  units  that  were  withdrawn.  The  upper  limit  is  for 

simplicity  and  introduces  no  technical  difficulty,  since  n(u)  =  0 
after  observation  ends.B 
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The  likelihood  (1.16)  applies  for  any  absolutely  continuous  life 
distribution.  In  the  important  special  case  of  an  exponential  life 
distribution  model,  f(x  |  X)  *  Xe  ^X  ,  the  likelihood  of  the  observed 
outcome  takes  the  simpler  form: 


,  k 

X  exp 


(1.17) 


L(X  I  D)  * 


<30 

-X  J  n 


(u)du 


exp 


oo 

-X  j"  n 


(u)du 


,  k  >  1 


k  -  0 


The  following  theorem  is  obvious  from  (1.17). 


1.8  Theorem 

Assume  that  the  test  plan  satisfies  Assumptions  (1),  (2)  and  (3)  of 
the  General  Sampling  Plan.  Assume  that  k  failures  and  the  number  of 
units  operating  at  age  u  ,  n(u)  ,  u  >_  0  ,  are  observed  and  that  the 
model  is  the  exponential  density  f(x  |  X)  =*  Xe  .  Then 


(a)  k  and  T  *  J  n(u)du  together  is  a  sufficient  statistic  for  X 

0 

(b)  I  *-s  the  MLE  for  X  . 

Note  that  the  MLE,  ^  >  f°r  *  represents  the  number  of  observed 
failures  divided  by  the  total  time  on  test.  This  conforms  with  the 
results  obtained  for  the  MLE  under  all  the  test  plans  consic-red  in 


Chapter  1. 
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The  maximum  likelihood  estimator  is  the  mode  of  the  posterior 
density  corresponding  to  a  uniform  prior  (over  an  interval  containing 
the  MLE) .  A  uniform  prior  is  often  a  convenient  reference  prior. 

Under  suitable  circumstances,  the  analyst's  actual  posterior  distri¬ 
bution  will  be  approximately  what  it  would  have  been  had  the  analyst's 
prior  been  uniform.  To  ignore  the  departure  from  uniformity,  it  is 
sufficient  that  the  analyst's  actual  prior  density  changes  gently  in 
the  region  favored  by  the  data  and  also  that  the  prior  density  not 
too  strongly  favors  some  other  region.  This  result  is  rigorously  ex¬ 
pressed  in  the  Principle  of  Stable  Estimation  (see  Edwarls,  Lindman 
and  Savage  (1963)).  DeGroot  (1970),  pages  198-201,  refers  to  this 
result  under  the  name  of  precise  measurement. 

1.9  Example 

The  exact  likelihood  can  be  calculated  explicitly  for  specified 
stopping  rules.  Suppose  that  withdrawal  times  are  determined  in 
advance.  Then  the  likelihood  is 


(1.18)  L(r(u)  ,  u  _>  0  |  D)  = 

where  n(xs)  is  ctle  number  surviving  just  prior  to  the  observed  failure 
at  age  xg  .  To  see  this  consider  the  airplane  engine  data  in  Example  1.4. 
Using  Figure  1.3  as  a  guide,  the  likelihood  will  have  the  following  factors 


k 

n 

S“1 


n(xs)r(xs) 


-J  n(u)r (u)du 
0 


1.  For  the  interval  [O.x^^j  we  have  the  contribution 
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4r(x^)  exp 


(1) 

-  J 


4r (u)du 


corresponding  to  the  probability  that  all  4  units  survive  to 
x^  and  the  first  failure  occurs  at  x^  » 

2.  For  the  interval  (x^j  ’^(1)  ^  we  ^ave  the  contribution 


exp 


'(1) 


-  J 


3r (u)du 


L(l) 


corresponding  to  the  conditional  probability  that  the  remaining 
3  survive  this  interval; 

3.  For  the  interval  (£^,x^]  we  have  the  contribution 


^r<X(2) >  exp 


(2) 


2r (u)du 


'Cl) 


corresponding  to  the  conditional  probability  that  the  remaining 
2  units  survive  to  x(2)  and  a  failure  occurs  at  x(2)  ; 

4.  For  the  interval  ^x(2)’^(2)^  we  ^ave  contribution 


exp 


(2) 

J*  r(u)du 


“(2) 


corresponding  to  the  conditional  probability  that  the  remaining 
unit  survives  to  age  i ^  .  Multiplying  together  these  condi¬ 
tional  probabilities,  we  obtain  a  likelihood  having  the  form 


shown  in  (1.18). 
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2.  PARAMETER  ESTIMATORS  AND  CREDIBLE  INTERVALS 

In  the  previous  section  we  saw  how  to  calculate  the  likelihood 
function  for  general  life  distributions.  This  is  required  in  order 
to  calculate  the  posterior  distribution.  Calculation  and  possibly 
graphical  display  of  the  posterior  density  would  conceivably  com¬ 
plete  our  data  analysis. 

If  we  assume  a  life  density  p(x  |  6)  and  tt (0 )  is  the  prior, 
then  p(x,0)  *  p(x  |  8)tt(9)  is  the  joint  density  and  p(x)  = 

J*  p(x  |  0)ir(0)d0  is  the  marginal  or  predictive  density.  Given  data 

D  and  the  posterior  density  ir(9  |  D)  ,  the  predictive  density  is 


p (x  j  D)  =  J  p(x  I  0 ) 7T (0  I  D)d8  . 

If  asked  to  give  the  probability  of  survival  until  time  t  ,  we  would 
calculate 


P(X  >  t  |  D)  = 


D)dx  . 


2 . 1  Example 

For  the  exponential  density  Xe  ,  k  observed  failures,  T 

total  time  on  test,  and  the  General  Sampling  Plan,  the  likelihood  is 
k  -XT 

proportional  to  X  e  .  For  the  natural  conjugate  prior, 


ir(X) 


,a,  a-1  -bX 
d  X  e 

r(a) 
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the  posterior  density  is 

it  (A  |  k  ,  T)  -  (b  +  T)a+kAa+k'1e"(b+T)X//r(a  +  k)  . 
In  this  case  the  probability  of  survival  until  time  t  is 


(2.1) 


P(X  >  t  k  ,  T) 


CO 


(A  I  k  ,  T)dA 


b  +  T  V 

b  +  t  +  T/ 


\a+k 


Bayes  Estimators 

We  will  need  the  following  notation: 


E[9] 


/ 


0TT(e)de 


and 


E[0  ]  D] 


07r  ( 9  |  D)d0  . 


Of  course,  E[0]  is  the  mean  of  the  prior  distribution  while  E [ 9  |  D] 
is  the  mean  of  the  posterior  distribution. 

We  wish  to  select  a  single  value  as  representing  our  "best" 
estimator  of  the  unknown  parameter  0  .  To  define  the  best  estimator 
we  must  specify  a  criterion  of  goodness  (or  equivalently,  of  poorness). 


23 


Statisticians,  being  pessimistic  by  nature,  measure  the  poorness  of 
an  estimator  9  by  the  expected  "loss"  resulting  from  their  estimator 
9  .  One  very  popular  loss  function  is  squared  error  loss:  specifically, 
having  observed  data  D  and  determined  the  posterior  density  ir(9  |  D)  , 
the  expected  squared  error  loss  i3  given  by: 

(2.2)  E[ (0  -  9) 2  |  D]  ; 

the  expectation  is  calculated  with  respect  to  the  posterior  density 
it (9  |  D)  .  We  choose  a  point  estimator  9  so  as  to  minimize  the  ex¬ 
pected  squared  error  loss  in  (2.2);  i.e.,  we  choose  9  to  satisfy: 

(2.3)  minimum  E[(9  -  a)2  |  D]  -  E[(9  -  9)2  |  D]  . 

a 

To  find  the  minimizing  value  9  ,  we  add  and  subtract  E(9  |  D) 
in  the  loss  function  to  obtain: 

E [ (9  -  a)2  |  D]  *  E[ (9  -  E(9  |  D))2  |  D]  +  [E(9  |  D)  -  a]2  . 

Since  we  wish  to  minimize  the  right  hand  side,  we  set  a  *  E(9  |  D)  , 
which  then  represents  the  solution  to  (2.3).  The  resulting  estimator, 

E(9  |  D)  ,  the  mean  of  the  posterior,  is  called  the  Bayes  estimator 
with  respect  to  squared  error  loss. 

2 . 2  Theorem 

The  Bayes  estimator  of  a  parameter  9  with  respect  to  squared 
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loss  is  the  mean  E(9  |  D)  of  the  posterior  density. 

Another  loss  function  in  popular  use  is  the  absolute  value  loss 
function: 

(2.4)  E[|e  -  e|  I  D]  . 

To  find  the  minimizing  estimator  using  this  criterion,  we  choose  0 
to  satisfy: 

(2.5)  minimum  E [ | 8  -  a|  j  D]  =  E[|0  -  0j  |  D]  . 

a 

It  is  easy  to  show: 

2.3  Theorem 

The  Bayes  estimator  of  a  parameter  0  with  respect  to  the  absolute 
value  loss  function  is  the  median  of  the  posterior  density.  Specifically, 

A 

the  estimator  0  satisfies: 

9  ® 

ir(0  |  D)d0  =  J  tt(0  |  D)d9  =  ~  . 

0 

The  proof  is  left  to  the  Exercises. 

Of  course,  the  prior  density  and  the  loss  function  enter  crucially 
in  determining  a  "best"  estimator.  However,  no  matter  what  criterion 
is  used,  all  the  information  concerning  the  unknown  parameter  0  is 
contained  in  the  posterior  density.  Thus,  a  graph  of  ir(0  j  D)  is  more 
informative  than  any  single  parameter  of  the  posterior  density,  whether 
it  be  the  mean,  the  median,  the  mode,  a  quartile,  etc. 


V 
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2.4  Example 

x 

1  *  Q 

Assume  that  lifetime  is  governed  by  the  exponential  model,  —  e 
Suppose  we  conjecture  that  E[6  |  k  ,  T]  ,  for  sampling  plan  with 
k  ,  T  sufficient,  is  linear  in  T  for  fixed  k  .  It  turns  out  that 
such  a  linear  relationship  holds  if  and  only  if  we  use  as  our  prior  the 
natural  conjugate  prior: 


r(6) 


_b 

bV(a+1)e~9 

r(a) 


(See  Section  3  of  Chapter  1  for  a  discussion  of  this  natural  conjugate 
prior  and  Diaconis  and  Ylvisaker  (1979)  for  a  proof  of  this  result  and 
for  more  general  results  of  this  kind.)  The  corresponding  Bayes 
estimator  with  respect  to  squared  error  loss  is: 


(2.7) 


E[9  |  k  ,  T] 


(b  +  T) 

(a  +  k  -  1) 


However,  the  natural  conjugate  prior  would  not  be  appropriate  if  we 
believed,  for  example,  that  9  could  assume  values  only  in  two  disjoint 
intervals.  Under  this  belief,  a  bimodal  prior  density  would  be  more 
natural,  and  the  corresponding  estimator  E [ 0  |  D]  would  very  likely 
be  difficult  to  obtain  in  closed  form  such  as  in  (2.7).  However 
E [ 0  |  D]  could  be  computed  by  numerical  integration. 

There  are  many  other  functions  of  unknown  parameters  for  which  we  may 
want  the  Bayes  estimator  with  respect  to  squared  error  loss.  For  example, 
we  may  wish  to  estimate  the  probability  of  survival  until  age  t  for  the 
exponential  model;  i.e.,  estimate 
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(2.8)  g(8)  -  exp  -  |-J  . 


It  is  easy  to  show  in  this  case  that 


(2.9) 

g  -  E 

exp 

t_ 

1  D 

0  . 

is  the  Bayes  estimator.  If  ir(0)  is  the  natural  conjugate  prior,  then 
it  is  easy  to  verify  that 

*  _  b  +  T  ja+k 
8  *  [b  +  t  +  t]  ’ 


i.e.,  this  is  the  Bayes  estimator  of  the  probability  of  survival  to  age 


given  total  time  on  test  T  and  k  observed  failures.  Note  that 
a+k 

is  precisely  the  marginal  probability  of  survival 


g 


b  +  T 


b  +  t  +  T 
until  time  t 


Credible  Intervals 

As  we  have  seen,  Bayes  estimators  correspond  to  certain  functions 
of  the  posterior  distribution  such  as  the  mean,  the  mode,  etc.  A 
credible  set  or  interval  is  another  way  of  presenting  a  partial  des¬ 
cription  of  the  posterior  distribution. 

Specifically,  we  choose  a  set  C  on  the  positive  axis  (since  we 
are  dealing  with  lifetime)  such  that 


J  tt(0  |  D)d0  -  1  -  a  . 

C 


(2.10) 
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Such  a  set  C  is  called  a  Bayesian  (1  -  a) 100  percent  credible  set 
(or  credible  interval  if  C  is  an  interval)  for  9  . 

Obviously,  the  set  C  is  not  uniquely  determined.  It  would  seem 
desirable  to  choose  the  set  C  to  be  as  small  (e.g.,  least  length, 
area,  volume)  as  possible.  To  achieve  this,  we  seek  a  constant 
c^_a  and  a  corresponding  set  C  such  that: 

(2.11)  C  =■  {9  |  ir(0  |  D)  >  ^  } 
and 

(2.12)  J  tt ( 0  I  D)d6  »  1  -  a  . 

C 

A  set  C  satisfying  (2.11)  and  (2.12)  is  called  a  highest  posterior 
density  credible  3et  [Box  and  Tiao  (1973) J.  In  general,  C  would 
have  to  be  determined  numerically  with  the  aid  of  a  computer. 

For  the  exponential  model  Xe  ,  we  have  seen  that  the  natural 
conjugate  prior  (Table  3.1,  Chapter  1)  is  the  gamma  density.  Since 
the  gamma  density  is  a  generalization  of  the  chi-square  density,  we 
recall  the  definition  of  the  latter  so  that  we  can  make  use  of  it  to 
determine  credible  intervals  for  the  failure  rate  of  the  exponential. 

2.5  Definition 

A  random  variable  >"(n)  having  density: 

x 
2 

—  for  x  >  0  ,  n  *  1,2,  ...  , 


(2.13) 


f  ,  (x) 

X  00 


H-  i 
2 

x  exp 


n 

22T 
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is  called  a  chi-square  random  variable  with  n  degrees  of  freedom  (d.f.). 

A  table  of  percentage  points  of  the  chi-square  distribution  may  be 
found  in  Pearson  and  Hartley  (1958).  In  addition,  chi-square  programs 
are  available  for  more  extensive  calculations  using  electronic  computers 
and  programmable  calculators. 

2 

It  is  easy  to  verify  that  the  x  random  variable  with  2n  d.f. 
is  distributed  as  2(Y1  +  Y£  +  —  +  YQ)  ,  where  Y^.Y^,  ...,  Yn  are 
independent,  exponentially  distributed  random  variables  with  mean  one. 
Thus,  we  obtain  the  following  result  useful  in  computing  credibility 
intervals  for  the  failure  rate  of  the  exponential  model  with  corresponding 
natural  conjugate  prior. 


2.6  Theorem 

Let  k  failures  and  total  time  on  test  T  be  observed  under 

sampling  as  suctions  (1),  (2)  and  (3),  (Section  1)  for  the  exponential 
“Ax  ** 

model  Ae  .  Let  A  have  posterior  density  corresponding  to  the 
natural  conjugate  prior 


it  (A) 


,a,a-l  -bA 
b  A  e 

r(a) 


with  a  an  integer.  Then 


(2.14) 


X  [2(a  +  k)  ]  X  [2 (a  +  k)  ] 

—  1  -  — 

2  :  x  2 

<  x  <  - - — — — - 


2(b  +  T)  - 


1  -  a  , 


2 

where  Xo(n)  is  the 


with  n  d.f.;  i. e. , 


1008  percentage  point  of  a  chi-square  distribution 

X2g(a) 

J  f  _  (x)dx  -  6  . 

0  x  (n) 
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Remark: 

2 

Because  of  the  lack  of  symmetry  of  the  x  density,  the  interval 
in  (2.14)  is  not  the  highest  posterior  density  credible  interval. 


Proof: 


It  is  easy  to  verify  that  (b  +  T)A  given  the  data  has  a  gamma  density, 


, a+k-1  -A 

A _ e 

r(a  +  k)  ’ 


corresponding  to  the  density  of  +  . . .  +  Ya+^  ,  where  the  Y's 
are  independent  unit  exponential  random  variables.  Hence 


2X(b  +  T)  -t  2(Y1+  ...  *Ya+k)  . 

where  ■  denotes  stochastic  equality;  i.e.,  2A(b  +  T)  has  a  chi- 
st 

square  density  with  2 (a  +  k)  d.f.  ■ 


2.7  Corollary 

For  2 (a  +  k)  large  (say  2(a  +  k)  >  30)  ,  the  normal  approximation 
provides  the  approximate  credibility  statement 


(2.15)  P 


(a  +  k)  +  (a  +  k)^z  (a  +  k)  +  (a  +  k)^z 


b  +  T 


<  A  < 


b  +  T 


1  D 

I  J 


1  -  a 


^  2 
where  z^  satisfies  J  4>(u)du  ■  a  and  4> (u)  ■  - 


u 

1  2  , 

e  is  the  normal 

SlU 


density  with  mean  0  and  variance  1. 
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Proof : 


Since  the 


X2(2n) 


random  variable  can  be  written  as 


X2(2n)  =  2(Y1  +  Y2  +  ...  +  Yq) 


where  Y.,Y„,  ...,  Y  are  independent  unit  exponentials,  the  Central 
1  2  n 

Limit  Theorem  [e.g.,  Hoel,  Port  and  Stone  (1971)]  applies.  Note  that 
Ex2(2n)  ■  2n  and  Var  [x2(2n)]  *  4n  .  Thus, 

x2(2n)  -  2n 

l/n 


is  approximately  normal  with  mean  0  and  variance  1  by  the  Central  Limit 
Theorem.  ■ 


2.8  Corollary 

Let  k  failures  and  T  total  time  on  test  be  observed  under  the 

General  Sampling  Plan  assumptions  (1) ,  (2)  and  (3)  (see  Section  1) , 

x 

1*0 

for  the  exponential  model  —  e  .  Let  0  have  the  natural  conjugate 
prior  with  integer  a  ,  then 


(2.16) 


2(b  +  T) 
a [2(a  +  k)] 

'7 


<  e 


<  2 (b  +  T) 

X2[2(a  +  k)] 

7 


1  -  a  . 


Proof; 

Since  9  has  the  natural  conjugate  prior  distribution  for  the  model 
_x 

-I  ~  9  '1 

9  e  ,  then  X  ■  —  has  the  natural  conjugate  prior  for  the  model 
9 
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Xe  .  (2.16)  follows  from  (2.14).  ■ 

2 . 9  Examp le 

We  use  the  data  in  Exercise  7,  Chapter  1,  to  illustrate  the 

calculation  of  the  credibility  limits  of  (2.16).  In  this  case,  n  »  13  , 

k  =  10  ,  and  T  =  275  weeks.  Figure  2.1  displays  prior  and  posterior  95% 

credibility  limits  for  selected  values  of  the  natural  conjugate  prior 

_x 

1  0 

parameters  for  the  exponential  model  —  e  .  For  the  improper  prior 

o 

ir (9)  *  -7  corresponding  to  a  *  b  *  0  we  cannot  calculate  a  95%  interval. 
However,  the  posterior  is  proper  in  this  case  and  we  can  calculate  a  95% 
credible  interval.  This  is  the  same  as  the  usual  sample  theory  95%  con¬ 
fidence  interval.  It  is,  however,  wider  than  posterior  credible  intervals 
corresponding  to  proper  priors.  The  credible  interval  corresponding  to 
the  prior  with  a  »  8  ,  b  *  120  was  computed  using  Corollary  2.7.  For 
this  relatively  small  number  of  failures  (k  »  10)  ,  the  choice  of  a 
prior  is  crucial  and  should  be  considered  rather  carefully. 

2.10  Example:  Credible  Limits  on  the  Failure  Rate  Average 

In  many  life  test  data  collections  only  the  numbers  of  failures  in 
specified  time  intervals  are  recorded,  and  not  the  more  informative  times 
of  failure.  For  example,  as  part  of  its  quality  assurance  program,  a 
company  producing  semiconductors,  routinely  selects  components  from  a 
product  batch  and  subjects  them  to  an  accelerated  life  test.  Speci¬ 
fically,  early  failure  is  induced  by  operating  the  components  in  a  high 
temperature  environment.  At  the  end  of  each  week,  components  are  examined 
and  those  which  do  not  meet  specifications  are  designated  as  failed. 


POSTERIOR 


57 


55  t 


43.76  T 


(a  *  0  ,  b  ■  0)  (a  =  4  ,  b  *  60)  (a  ■  8  ,  b  =  120) 


FIGURE  2.1 


COMPARISON  OF  95%  CREDIBLE  INTERVALS  FOR  9  (k  -  10  ,  T  -  275  weeks) 


PRIOR 
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In  Table  2.1,  the  results  of  testing  n  *  840  semiconductors  are  recorded 
for  a  3  week  period. 


TABLE  2.1 

SEMICONDUCTOR  LIFE  TEST  DATA 


Time  Interval,  hours 

Observed  Number  of  Failures 

0  -  168 

9 

168  -  336 

6 

336  -  504 

13 

ds  £ 

Let  N(t)  be  the  number  of  failures  in  [ 0 , t ]  and  p  =  F(t)  , 
the  probability  of  failure,  in  the  interval.  Then 


(2.17) 


P[N(t)  -  k  I  P]  - 


n  k,  . n-k 
k  P  (1  -  P) 


From  Table  2.1  it  is  evident  that  failures  are,  relatively  speaking, 
rare  events.  Thus,  a  Poisson  distribution  approximation  to  the  binomial 

may  be  justified;  i.e.,  P[N(t)  ■  k  j  p]  *  e  np  .  Note  also  that 

the  MLE  for  p  in  the  first  interval  [0,168]  is  p  =  0.01  .  This 
suggests  the  approximation 


p  =  1  -  exp 


/ 


r(u) du 


I 


r(u)du  . 


Let  X  ■  —  /  r(u)du  ,  the  failure  rate  averaae  in  the  first  time  interval. 
C  0 

(If  the  failure  rate  is  constant,  then  of  course  so  also  is  the  failure 
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rate  average.)  Using  these  approximations  we  obtain  for  the  interval 
<0,t]  ; 


L(X  |  k  ,  t) 


(n.X  t )  k  -nAt 
k!  6 


The  natural  conjugate  prior  is  easily  seen  to  be: 


it  (A) 


baAa_1  -bA 

r(a)  e 


It  follows  that  the  posterior  is  given  by: 


"(A  |  k  ,  t) 


(b  +  nt)a'l~kAa't'k~1  -(b+nt)  A 
r(a  +  k)  e 


For  example,  suppose  we  choose  as  our  hyper-parameters  for  the  prior 
distribution  a  «  4  ,  and  *  10  ^  or  b  *  4  *  10^  .  Using  Theorem  2.6, 
we  can  compute  95%  credible  limits  on  the  failure  rate  average.  For  the 
first  time  interval  of  168  hours,  we  have 


P  js .  70  x  10' 


X  025^a  + 
2(nt  +  b) 


<  A 


x_975[2(a  +  k)] 


1.44  x  10 


-4 


-  1  -  2 (nt  +  b) 

) 


J 


0.95 


If  we  judge  failure  rate  parameters  for  distinct  time  intervals  to  be 
independent,  then  it  follows  that  they  are  also  independent  with  respect 
to  the  joint  posterior  distribution.  Using  this  fact  a  joint  posterior 
distribution  for  failure  rates  for  all  three  time  intervals  can  be  easily 


calculated. 
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WEIBULL  DISTRIBUTION 


Whenever  possible,  the  choice  of  a  life  distribution  model  should 
be  based  on  the  underlying  failure  mechanisms.  Simple  structures  com¬ 
posed  of  statistically  independent  components  have  been  used  to  derive 
life  distribution  models  valid  when  the  number  of  structural  components 
is  very  large. 

Suppose  a  structure  of  n  components  fails  as  soon  as  k  com¬ 
ponents  fail.  For  example,  each  of  the  .'trands  whose  lifetimes  under 
stress  are  recorded  in  Table  1.1,  Chapter  1,  is  composed  of  approximately 
270  filaments  bound  together  by  an  epoxy.  If  strand  failure  coincides 

with  the  kth  failure  of  filaments,  then  structure  lifetime  is  the 

k,n 

kth  smallest  of  n  lifetimes.  If  also  component  lifetimes  are  judged 

identically  distributed  and  independent,  then  there  are  only  two  possible 

limiting  structure  life  distributions  in  the  sense  that  there  exist 

00  00 

sequences  of  normalizing  constants  {a  }  ,  ,  {X  }  ,  such  that  for  all 

n  n=i  n  n=l 

real  x  : 


lim  P{X  U  n  -  a  )  <  x} 
n-K=  n  k,n  n  - 


exists.  The  limit  is  either 


(3.1) 


U(x-a)r 


i  r  -uk-i,  ... 

(fc  T  1)~;  I  e  u  du,a,X>0,x>a>_0, 


exp  [X  (x- a)  ] 


(k  -  1)!  J 

(3.2)  0 


-u  k-1 

e  u  du,-00<x<“> 


-oo<a<oo,X>0 


36 


[Smirnov  (1952)].  In  both  cases  a  is  a  location  parameter  and  X  is 
a  scale  parameter  while  a  and  k  are  shape  parameters. 

If  k  ■  1  ,  then  (3.1)  becomes 

(3.1')  W(x  |  a  ,  X  ,  a)  ■  1  -  exp  { —  [ X (x  -  a)]a}  ,  x  >_  a  >_  0  , 

the  Weibull  distribution ,  and  (3.2)  becomes 

(3.2')  A(x  )  a  ,  X)  ■  1  -  exp  {-e^X  }  ,  -®  <  x  <  00  . 

Thus,  if  X  is  the  structure  lifetime,  then  either  X  or  exp  (X)  has 
a  Weibull  distribution.  The  failure  rate  for  the  Weibull  distribution 
of  (3.1*)  is 


ry(x)  ■  atXa(x  -  a)0  1  for  x  >_  a  , 

and  0  elsewhere.  In  the  second  case  it  is 
(3.4)  r^(x)  *  *  exp  ~  a)  J  • 

For  all  parameter  values,  (3.4)  is  increasing  in  x  .  Hence,  if  we  wish 
to  allow  the  possibility  that  the  failure  rate  mav  be  decreasing  we  must 
choose  the  Weibull  model,  (3.1'),  with  a  <  1  . 

The  Weibull  model  appears  to  furnish  an  adequate  fit  for  some  strand 
lifetime  data  with  estimated  values  of  a  less  than  4.  On  the  other 
hand,  it  has  been  empirically  observed  that  for  strength  data,  estimates 
for  a  using  the  Weibull  model  are  often  large  (>  27  in  some  cases). 
This  suggests  that  (3.2')  may  provide  a  better  model  for  strand  strength 


data. 
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Inference  for  the  Weibull  Distribution 

The  Weibull  life  distribution  model  has  three  parameters:  a  ,  X  , 
and  a  .  The  parameter  a  >  0  is  a  threshold  value  for  lifetime;  before 
time  a  we  expect  to  see  no  failures.  If  there  is  no  physical  reason 
to  justify  a  positive  threshold  value,  the  analyst  should  use  the  two 
parameter  Weibull  model.  The  most  simple  model  compatible  with  prior 
knowledge  concerning  physical  processes  will  often  provide  the  most  in¬ 
sight.  The  Weibull  density  is 


c\  cr  l  i  \  ,<*,  \a-l  -[X(x-a)]a 

(3.5)  f(x  |  a  ,  a  ,  X)  ■  aX  (x  -  a)  e 

for  x  ^  a  and  0  elsewhere. 

Usually  we  wish  to  quantify  our  uncertainty  about  a  particular 
aspect  of  the  life  distribution,  such  as  the  probability  of  surviving 
x  hours.  For  the  three  parameter  Weibull  model,  this  is  given  by: 


(3.6) 


F(x  |  a  ,  X  ,  a)  *  exp  {-[X(x  -  a)]a}  . 


It  is  clearly  sufficient  to  assess  our  uncertainty  concerning  a  ,  X  , 
and  a  . 

Suppose  data  are  obtained  under  the  General  Sampling  Plan  (Section  1) . 
Let  x^x^,  ...,  x^  denote  the  unordered  observed  failure  ages  and  n(u) 
the  number  surviving  until  age  u  .  Then  by  Theorem  1.6  in  Section  1, 
the  likelihood  is  given  by: 


L(a  ,  a  ,  X  |  D) 


k,  ka 
a  X 


n  (x  -  a) 
i-1 


a-1 


-X 


ao 

I 


an(u) (u  -  a)a  ^du 


JJ 


exp 
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for  a  £  and  a  ,  X  >  0  .  Suppose  there  are  m  withdrawals  and  we 
pool  observed  failure  and  loss  times  and  relabel  them  as: 


0  '  C(0)  -  C(l)  -  C(2)  - 


i  t(k+m)  i  C 


Then,  for  a  ^  x.^  ,  i  -  1,2,  . ..,  k  ,  we  have 


QD 

I" 


(u) (u  -  a)a  1du 


C(i) 


k+m  y' 

l  (n  -  i  +  1)  I  (u  -  a)a 

i-1  J 


du  +  (n  -  k  -  m) 


L 

J  (u  -  a)a  ^du 


'  (k-Hn) 


Observation  is  confined  to  the  age  interval  [0,t]  . 

Two  important  deductions  can  be  made  from  (3.7): 

1.  The  only  sufficient  statistic  for  all  three  parameters  (or 
for  a  and  A  alone  when  a  ■  0)  is  the  entire  data  set. 

2.  No  natural  conjugate  family  of  priors  is  available  for  all 
three  parameters  (or  for  a  and  X  alone  when  a  »  0) . 
Consequently,  the  posterior  distribution  must  be  computed 
using  numerical  integration.  [see  Diaconis  and  Ylvisaker 
(1979)]. 

For  most  statistical  investigations,  a  and  perhaps  also  a 
would  be  considered  nuisance  parameters.  By  matching  our  joint  prior 
density  on  a  ,  X  and  a  with  the  likelihood  (3.7),  we  can  calculate 
the  posterior  density,  it  (a  ,  X  ,  a  |  D)  .  For  example,  if  a  is 
considered  a  nuisance  parameter,  then  we  would  calculate  the  marginal 
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density  on  A  and  a  as: 


A  |  D)da  . 


Credibility  Regions  for  Two  Parameter  Models 


Let  n(a  ,  A  |  D)  be  the  posterior  density  for  a  two  parameter 
model  such  as  the  Weibull  model  above  with  scale  parameter  A  and 
shape  parameter  a  .  To  find  the  so-called  "highest  posterior  density" 
credibility  region  for  a  and  A  simultaneously  (Section  2),  we  find 
a  constant  c(8)  by  sequential  search  such  that: 


(3.8)  R  -  [(a  ,  A)  |  ir(a  ,  A  |  D)  >  c(6) ] 


and 


it  (a  ,  A  j  D)dadA 


8  . 


The  region  R  defined  above  is  a  8(100)  percent  credibility  region 
for  a  and  A  .  For  unimodal  densities  such  regions  are  bounded  by 
a  single  closed  curve  C  which  does  not  intersect  itself  (i.e.  ,  a 
"simply  connected  region"). 

To  illustrate  the  use  of  Weibull  credibility  regions  we  have  computed 
credibility  regions  corresponding  to  the  data  in  Tables  3.1  and  3.2. 
Twenty-one  pressure  vessels  were  put  on  life  test  at  68%  of  their  ultimate 
mean  burst  stress.  A  pressure  vessel  is  filled  with  a  gas  or  liquid  and 
provides  a  source  of  mechanical  energy.  They  are  used  on  space  satellites 
and  other  space  vehicles.  After  13,488  hours  of  testing,  5  failures  were 


V 
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TABLE  3.1 

ORDERED  FAILURE  AGES  OF  PRESSURE 
VESSELS  LIFE  TESTED  AT  68%  OF  MEAN  RUPTURE 
STRENGTH  (n  -  21  ,  OBSERVATION  TO  13,488  HOURS). 


Number  of 
Failure 

Age  at  Failure 
(hours) 

1 

4000 

2 

5376 

3 

7320 

4 

8616 

5 

9120 

TABLE  3.2 

ORDERED  FAILURE  AGES  OF  PRESSURE  VESSELS 
LIFE  TESTED  AT  68%  OF  MEAN  RUPTURE  STRENGTH 
(FAILURES  BETWEEN  13,488  HOURS  AND  20,568  HOURS) 


Number  of 

Age  at  Failure 

Failure 

(hours) 

1 

14400 

2 

16104 

3 

20231 

4 

20233 
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recorded.  After  an  additional  7080  hours  of  testing,  an  additional  4 
failures  were  recorded. 

Figure  3.1  displays  credibility  contours  for  a  and  X  after 
13,488  hours  of  testing  and  again  after  20,568  hours  of  testing.  The 
posterior  densities  were  computed  relative  to  uniform  priors.  The 
posterior  density  computed  after  20,568  hours  could  also  be  inter¬ 
preted  as  the  result  of  using  the  posterior  (calculated  on  the  basis  of 
Table  3.1  and  a  flat  prior)  as  the  new  prior  for  the  data  in  Table  3.2. 
A  qualitative  measure  of  the  information  gained  by  an  additional  year 
of  testing  can  be  deduced  by  comparing  the  initial  (dark)  contours 
and  the  tighter  (light)  contours  in  Figure  3.1. 

To  predict  pressure  vessel  life  at  the  68%  stress  level,  we  can 
numt  rically  compute 


P[X  >  t  |  D] 


ff 


e-<*t) 


a 

tt(o  ,  X  |  D)dadX 


where  ir(a  ,  X  j  D)  must  be  numerically  computed  using  the  given  data, 
D  . 

If  the  mean  life 


9 


or  the  standard  deviation  of  life  are  of  interest,  their  posterior 
densities  can  be  computed  by  making  a  change  of  variable  and  integrating 
out  the  nuisance  parameter.  For  example,  if  a  ■  0  in  the  Weibull 
model  and  we  are  interested  in  the  mean  life,  9  ,  we  can  use  the 
Weibull  density  in  terms  of  a  and  6 


(per  1000  days) 
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0.40  0-00  1.20  1.60  2.00  2.40  2.00  9.20  9.60  4.00  4.40  4.00 


FIGURE  3.1 

HIGHEST  PROBABILITY  DENSITY  CONTOURS  FOR  a  AND  1  FOR 
KEVLAR/EPOXY  PRESSURE  VESSEL  LIFE  TEST  DATA.  THE  PRESSURE 
VESSELS  WERE  TESTED  AT  68Z  STRESS  LEVEL. 
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to  compute  the  joint  posterior  density  it  (a  ,  0  |  D)  .  The  prior  for 
a  and  X  must  be  replaced  by  the  induced  prior  for  a  and  0  . 

This  may  be  accomplished  by  a  change  of  variable  and  by  computing  the 
appropriate  Jacobian.  The  marginal  posterior  density  of  0  is  then 

GO 

ir(0  |  D)  *  J*  it  (a  ,  0  j  D)da  . 

0 

This  can  then  be  used  to  obtain  credibility  intervals  on  0  . 
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4.  EXERCISES 

1.  Twenty-one  pressure  vessels  are  subjected  to  a  static  stress 
equivalent  to  68%  of  the  mean  rupture  stress.  Testing  begins 
June  22,  1977.  Failures  as  of  November  1,  1979  are  recorded 
in  Table  4.1  in  hours. 


TABLE  4.1 

PRESSURE  VESSEL  FAILURE  TIMES  IN  HOURS 

4000 

5376 

7320 

8616 

9120 

14400 

16104 


Use  a  time  transformed  exponential  model  with  a  »  1.3  (see 
Example  1.1).  Plot  the  posterior  density  of  X3  relative 
to  a  uniform  prior.  Let  Rq(x)  “  ,  what  is  the  MLE? 

The  mean  life  for  this  model  is 


9  -  r (1  +  l/a)/X1/a  . 


What  is  the  mode  for  the  corresponding  posterior  density  for  9  ? 

2.  The  probability  densities  (1.12)  and  (1.13)  were  calculated 

for  the  case  where  unit  i  is  observed  until  either  it  fails 

at  age  x^  or  is  withdrawn  from  observation  at  age  .  Suppose 

we  have  only  ordered  observations  x,,.  <  x,„.  <  ...  <  x,, . 

(1)  -  (2)  -  -  (k) 

and  <  ...  <  Z.  .  ;  i.e.,  we  no  longer  know  which 

(1)  -  (2)  -  -  (m)  6 

unit  fails  at  age  or  Was  withdrawn  from  observation  at 
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age  •  Recalculate  the  likelihood  (1.12)  ana  (1.13)  for 

this  case.  [Hint:  see  Example  1.9.] 

3.  Suppose  only  the  ordered  failure  ages  1.  x(2)  —  *  *  *  —  x(k) 

and  the  ordered  ages  at  withdrawal  —^{2 )  —  —  ^(m) 

are  given.  Furthermore,  suppose  withdrawal  age  W  is  random 
with  withdrawal  rate  p(*)  ,  the  analogue  of  failure  rate. 

Show  that  the  exact  expression  for  the  likelihood  (1.16)  is 
now 


L(r(u)  ,  p(u)  ,  u  >  0  I  D) 


m 


CO 

-  J  n(u)p (u)du 


EV 

‘  ^  "Ks>)r(,W  *** 


00 

-  J"  n(u)r(u)du 


where  n(u  )  is  the  number  surviving  just  prior  to  age  u 
4.  Let  the  failure  rate  function  be 


r(x) 


\\1  ,  0  <_  x  <.  tx 


fx2  ’  x  -  cl  ' 


Given  n(u)  ,  the  number  surviving  until  time  u  >_  0  and 
observed  failure  times 


x(l)  -  x(2)  -  *•'  -  x(k)  * 

use  (1.16)  to  compute  the  likelihood  up  to  a  constant  factor 
independent  of  and  X2  •  Assume  noninfo rmative  stopping 
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rules  are  in  effect  and  that  t^  is  known. 

5.  Define  the  pth  percentile,  xp  ,  as  the  solution  to 


r 

J  f(x  I  8)dx  ■  p  . 

0 

For  the  exponential  model  f(x  j  9)  ■  j  exp  and  the 

natural  conjugate  prior  on  9  ,  find  the  Bayes  estimator 
for  Xp  corresponding  to  squared  error  loss.  Assume  that 
k  failures  and  T  total  time  on  test  have  been  observed. 

If,  for  example,  p  «  0.01  then  x^  might  qualify  as  the 
warranty  period.  The  warranty  would  then  be  valid  for 
x  hours  (months,  years,  etc.) 

6.  Let  Xp(9)  be  the  pth  percentile  for  the  exponential  model 
f(x  |  9)  »  exp  "f )  as  described  in  Exercise  5.  For  the 
data  in  Exercise  7,  Chapter  1,  compute  an  upper  99%  credible 
interval  on  x^  .  Use  the  same  natural  conjugate  prior 
specified  in  Exercise  7,  Chapter  1. 

7.  Show  that  minimum  E  [ 1 9  -  a  I ]  is  attained  when  a  satisfies 

<1  TT  1 


CL 

j*  ir(6)d6  -  j  ; 


i.e.,  a  is  the  median  for  tt  .  [The  same  result  holds  for 
the  posterior  density.  See  (2.6)). 

8.  Using  the  data  in  Table  2.1,  calculate  95%  credibility  intervals 
on 
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168 

X1  ‘  168  /  r<")d"  • 

0 

336 

X2  ■  168  J  rCu)d“  ' 
168 


and 


504 

X3  =  158  /  r(u)du  * 

336 

3 

Use  natural  conjugate  priors  with  a  ■  4  and  b  =  4  *  10 
for  both  A^  ,  A^  and  A^  .  Suppose  we  judge  them  indepen¬ 
dent.  Initially,  840  semiconductors  were  put  on  life  test. 
10.  Let  F (x  |  a  ,  A)  *  1  -  exp  [-(Ax)3]  .  Assume  a  is  known 
and  let  y  =  Aa  .  Let 


*(y) 


,  a  a-1  -by 
bye 

T(a) 


Calculate  the  Bayes  estimator  with  respect  to  squared  error  loss 
of  the  probability  of  survival  until  age  x  given  k  failures 
and  n(u)  ,  n  >_  0  where  n(u)  is  the  number  of  units  surviving 
until  age  u  . 

-Ax 

11.  Assume  the  exponential  model  Ae  and  suppose  k  failures 
and  T  total  time  on  test  have  been  observed.  Using  the 

baAa-1e"bX 

natural  conjugate  prior  tt(A)  *  - F(a) -  *  coniPute  the  posterior. 

Now  suppose  that  the  resulting  posterior  is  used  as  a  prior  with 
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the  same  data,  i.e.,  no  additional  data  are  collected  to  compute 
a  new  posterior.  Furthermore,  suppose  that  this  process  is 
repeated  infinitely  often  with  the  same  data. 

(a)  What  wi±l  be  our  final  conclusion? 

(b)  What  is  the  fallacy  in  this  inference  procedure? 

12.  Suppose  n  units  are  put  on  life  test.  The  stopping  rule  is 
as  follows.  If  the  first  failure  occurs  before  time  t^, 
stop  at  the  time  of  first  failure.  Otherwise  continue  testing 
until  time  t£  >  t  .  Compute  the  likelihood  for  the  exponen¬ 
tial  model.  Is  the  stopping  rule  noninformative?  Why? 

13.  Consider  the  following  life  test  for  n  +  1  units,  each  with 
life  distribution  F(*  j  8)  .  A  unit  is  selected  at  random 
and  life  tested  to  failure,  say  at  time  x  .  This  failure 
time  is  then  used  to  provide  the  stopping  rule:  test  n  units 
to  time  x  .  Is  this  stopping  rule  informative  relative  to  8  ? 


5.  NOTES  AND  REFERENCES 
Section  1 

In  the  General  Sampling  Plan  we  needed  to  assume  that  any  stopping 
rules  used  were  noninformative  concerning  the  failure  distribution.  The 
need  for  this  assumption  was  pointed  out  by  Raiffa  and  Schlaiffer  (1961). 
Examples  of  informative  stopping  rules  were  given  by  Roberts  (1967)  in 
the  context  of  two  stage  sampling  of  biological  populations  to  estimate 
population  size  (so-called  capture-recapture  sampling) . 

Section  2:  Unbiasedness 

The  posterior  mean  is  a  Bayes  estimator  of  a  parameter,  say  0  , 
with  respect  to  squared  error  loss.  It  is  also  a  function  of  the  data. 

An  estimator,  0(D)  ,  is  called  unbiased  in  the  sample  theory  sense  if 

Ep[0 (D)  |  0]  -  B 

for  each  0  E  0  .  No  Bayes  estimator  (based  on  a  corresponding  proper 
prior)  can  be  unbiased  in  the  sample  theory  sense  [Bickel  and  Blackwell 
(1967)]. 

Most  unbiased  estimators  are  in  fact  inadmissible  in  the  sample 

theory  sense  with  respect  to  squared  error  loss.  For  example,  0(D)  =* 

is  a  sample  theory  unbiased  estimator  for  the  mean  of  the  density 
_x 

1  ^  0 

—  e  under  sample  plan  (a)  of  Chapter  1.  However  it  is  inadmissible 
in  the  sense  that  there  exists  another  c0(D)  with  c  f  1  such  that, 
for  all  0 


Ep[ [cQ (D)  -  0]2  [  0]  <  Ep[[0(D)  -  0]2  |  0]  . 
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To  find  this  c  ,  consider  Y  *  and  note  EY  *  1  •  Then  we 

need  only  find  c  such  that 

Ep[(cY  -  l)2  |  9] 

EY 

is  minimum.  This  occurs  for  c„  =  — r  which  is  clearly  not  1.  Hence 

U  EY^ 

0(D)  is  sample  theory  inadmissible.  Sample  theory  unbiasedness  is 
not  a  viable  criterion. 

For  large  k  ,  0 (D)  m  —  will  be  approximately  the  same  as  our  Bayes 
T 

estimator.  However,  ^  is  not  recommended  for  small  k  . 

Since  tables  of  the  chi-square  distribution  have  in  the  past  been 
more  accessible  than  tables  of  the  gamma  distribution,  we  have  given 
the  chi-square  special  treatment.  However  with  modern  computing  facili¬ 
ties,  we  really  only  need  to  use  the  more  general  gamma  distribution. 

Confidence  Intervals 

A  (1  -  o)100%  confidence  interval  in  the  sample  theory  sense  is  one 
such  that  if  the  experiment  is  repeated  infinitely  often  (and  the  interval 
recomputed  each  time)  then  (1  -  a) 100%  of  the  time  the  interval  will 
cover  the  fixed  unknown  true  parameter  9  .  Since  confidence  intervals 
do  not  produce  a  probability  distribution  on  the  parameter  space  for  0  , 
they  cannot  provide  the  basis  for  action  in  the  decision  theory  sense; 
i.e.,  a  decision  maker  cannot  use  a  sample  theory  confidence  interval  to 
compute  an  expected  utility  function  which  can  then  be  maximized  over  his 

set  of  possible  decisions. 

-Xx  X 

If  for  Xe  we  choose  the  improper  prior,  ir(X)  -  —  ,  then  the 

A 

chi-square  (1  -  a) 100%  credible  intervals  and  the  sample  theory 


l 
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(1  -  a) 100%  confidence  intervals  agree.  Unfortunately,  such  improper 
credible  intervals  can  be  shown  to  violate  certain  rules  of  logical 
behavior.  Lindley  (personal  communication)  provides  the  following 
simple  illustration  of  this  fact  for  the  exponential  model  Xe 
Suppose  n  units  are  put  on  test  and  we  stop  at  the  first  failure,  so 
that  T  =  ^(i)  •  ^ow  ^  8*ven  ^  also  has  density  Xe  so  that 


(In  2) 
T 


is  a  50%  improper  upper  credible  limit  on 


X 


i. e. , 


(5.1) 


X  < 


(In  2) 


-  T 


it(X) 


1 

X 


0.50  . 


Suppose  now  that  T  is  observed  and  we  accept  the  probability 
statement  (5.1).  Consider  the  following  hypothetical  bet. 


(i) 

If 

X  < 

(In  2) 

T 

-T 

we  lose  the  amount  e  ; 

(ii) 

If 

X  >_ 

(In  2) 

T 

-T 

we  win  e 

We  can  then  pretend  that  the  true  X  is  somehow  revealed  and 
bets  are  paid  off.  If  we  believe  statement  (5.1),  then  given  T 
such  a  bet  is  certainly  fair. 

Now  let  us  compute  our  expected  gain  before  T  is  observed 
(preposterior  analysis).  This  is  easily  seen  to  be  (conditional  on  X) 


In  2 

X 


I 


Xe  ^Ce  Cdt  +  J*  Xe  ^Ce  Cdt 

In  2 

X 


X 

1  +  X 


1 


which  is  negative  for  all  X  >  0  .  Note  that  this  is  what  we  subjectively 
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expeat,  since  as  (improper)  Bayesians,  every  probability  (and  presum¬ 
ably  even  an  improper  prior)  is  subjective. 

The  contradiction  lies  in  the  observation  *-hat 

1.  Conditional  on  X  and  prior  to  observing  T  ,  our  expected 
winnings  are  negative  for  all  X  ; 

2.  Conditional  on  T  ,  our  expected  loss  is  zero  (using  the 
improper  prior  i»(X)  »  y). 

The  source  of  the  contradiction  is  that  we  have  not  measured  our 

uncertainty  for  all  events  by  probability.  For  example,  we  have 

assigned  the  value  *  to  the  event  A  <  Xq  for  all  \Q  >  0  ;  i.e., 

X  X 

o  o 

/  ir(X)dX  -  /  y  dA  »  ®  .  See  Appendix  B  for  a  proof  that  for  any  set 
0  0 

of  uncertainty  statements  that  are  not  probabilistically  based  (relative 
to  proper  distributions),  a  system  of  bets  can  be  constructed  which  will 
result  in  the  certain  loss  of  money.  A  bet  consists  of  paying  pz  <  z 
dollars  to  participate  with  the  understanding  that  if  an  event  E  occurs 
you  win  z  dollars  and  otherwise  you  win  nothing. 

Section  3 

The  Weibull  distribution  is  one  of  several  extreme  value  distri¬ 
butions.  See  Barlow  and  Proschan  (1975),  Chapter  8,  for  a  more  advanced 
discussion  of  extreme  value  distributions. 
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